Modularbeit Wiener Filter¶

Herleitung des Wiener Filters:¶

Image-Restauration

Suche: $$\begin{equation*} \underset{W}{argmin} \; \mathbb{E}\Big[\big\|\hat f - f\big\|^2\Big] \end{equation*}$$

Annahmen:¶

$$\begin{align*} g = h \ast f + n \; &\Longrightarrow \; G = H \cdot F + N \tag{1}\\ \hat f = g \ast w \; &\Longrightarrow \; \hat F = G \cdot W \tag{2}\\ \end{align*}$$

$$\begin{align*} S_{FF} &= \mathbb{E}\Big[\big\|F\big\|^2\Big] \tag{3}\\ S_{NN} &= \mathbb{E}\Big[\big\|N\big\|^2\Big] \tag{4}\\ \end{align*}$$

Das Bild und das Rauschen sind unkorreliert: $$\begin{align*} S_{FN} = \mathbb{E}\Big[F^*N\Big] = \mathbb{E}\Big[FN^*\Big] = 0 \tag{5}\\ \end{align*}$$

Herleitung:¶

Minimierung des quadratischen Fehlers vom Orginalbild und des rekonstruierten Bildes: $$\begin{equation*} \mathbb{E}\Big[\big\|\hat F - F\big\|^2\Big] \overset{(2)}{=} E\Big[\big\|GW - F\big\|^2\Big] \\ \end{equation*}$$

$$\begin{equation*} = \mathbb{E}\Big[\big\|GW\big\|^2 + \big\|F\big\|^2 - (GW)^*F - GWF^*\Big] \\ \end{equation*}$$

$$\begin{equation*} = \big\|W\big\|^2 \,\mathbb{E}\Big[\big\|G\big\|^2\Big] + \mathbb{E}\Big[\big\|F\big\|^2\Big] - W^* \,\mathbb{E}\Big[G^*F\Big] - W \,\mathbb{E}\Big[GF^*\Big] \\ \end{equation*}$$

$$\begin{equation*} \overset{(1,3)}{=} \big\|W\big\|^2 \,\mathbb{E}\Big[\big\|HF + N\big\|^2\Big] + S_{FF} - W^* \,\mathbb{E}\Big[(HF + N)^*F\Big] - W \,\mathbb{E}\Big[(HF + N)F^*\Big] \\ \end{equation*}$$

$$\begin{equation*} = \big\|W\big\|^2 \,\mathbb{E}\Big[\big\|HF\big\|^2 + \big\|N\big\|^2 + (HF)^*N + HFN^*\Big] + S_{FF} - W^* \,\mathbb{E}\Big[H^*\big\|F\big\|^2 + N^*F\Big] - W \,\mathbb{E}\Big[H\big\|F\big\|^2 + NF^*\Big] \end{equation*}$$

$$\begin{align*} = &\big\|W\big\|^2 \Big(\big\|H\big\|^2 \,\mathbb{E}\Big[\big\|F\big\|^2\Big] + \mathbb{E}\Big[\big\|N\big\|^2\Big] + H^* \,\mathbb{E}\Big[F^*N\Big] + H \,\mathbb{E}\Big[FN^*\Big]\Big) + S_{FF} \\ &- W^* \Big(H^* \,\mathbb{E}\Big[\|F\|^2\Big] + \mathbb{E}\Big[N^*F\Big]\Big) - W \Big(H \,\mathbb{E}\Big[\big\|F\big\|^2\Big] + \mathbb{E}\Big[NF^*\Big]\Big) \\ \end{align*}$$

$$\begin{equation*} \overset{(3,4,5)}{=} \big\|W\big\|^2 \Big(\big\|H\big\|^2 S_{FF} + S_{NN} + H^*0 + H0\Big) + S_{FF} - W^* \Big(H^* S_{FF} + 0\Big) - W \Big(H S_{FF} + 0\Big) \\ \end{equation*}$$

$$\begin{equation*} = \big\|W\big\|^2 \Big(\big\|H\big\|^2 S_{FF} + S_{NN}\Big) + S_{FF} - W^*H^*S_{FF} - WHS_{FF} \end{equation*}$$

Ableiten und Nullstellen: $$\begin{align*} \frac{\partial \, \mathbb{E}\Big[\big\|\hat F - F\big\|^2\Big]}{\partial \, W} &= \frac{\partial \, \big\|W\big\|^2 \Big(\big\|H\big\|^2 S_{FF} + S_{NN}\Big) + S_{FF} - W^*H^*S_{FF} - WHS_{FF}}{\partial \, W} \\ &= W^* \Big(\big\|H\big\|^2 S_{FF} + S_{NN}\Big) - HS_{FF} \overset{!}{=} 0 \\ \end{align*}$$

$$\begin{align*} &\Longleftrightarrow W^* \Big(\big\|H\big\|^2 S_{FF} + S_{NN}\Big) = HS_{FF} \\ &\Longleftrightarrow W \Big(\big\|H\big\|^2 S_{FF} + S_{NN}\Big) = H^*S_{FF} \\ &\Longleftrightarrow W = \frac{H^*S_{FF}}{\|H\|^2 S_{FF} + S_{NN}} = \frac{H^*}{\|H\|^2 + \frac{S_{NN}}{S_{FF}}} \end{align*}$$

Implementation:¶

In [ ]:
from PIL import Image
import numpy as np
from matplotlib import pyplot as plt
from tqdm import tqdm
import math
In [ ]:
def read_image(image_path):
    image = np.array(Image.open(image_path).convert('L'))
    normalized_image = image / image.max()
    return normalized_image

def show_image(image):
    plt.imshow(image, cmap='gray')
    plt.axis('off')
    plt.show()
In [ ]:
def fft(image):
    return np.fft.fft2(image)

def ifft(image_fq):
    return np.abs(np.fft.ifft2(image_fq))

Mean Squared Error:¶

$$\begin{equation*} MSE = \mathbb{E}\Big[\big\| \hat f - f \big\|^2\Big] \end{equation*}$$

Peak Signal To Noise Ratio:¶

$$\begin{equation*} PSNR = 10\log_{10}\Bigg(\frac{MAX_f^2}{MSE}\Bigg) \end{equation*}$$

Accuracy Amoothness Error:¶

$$\begin{equation*} ASE = \mathbb{E}\Big[\big\| \hat g - h \ast \hat f \big\|^2\Big] + \lambda \, \mathbb{E}\Big[\big\| \hat f^{\prime} \big\|^2\Big] \end{equation*}$$

In [ ]:
def mse(f, f_hat):
    return (np.square(f - f_hat)).mean()

def psnr(f, f_hat):
    mse_value = mse(f, f_hat)
    if mse_value == 0:
        return 0
    return 10 * np.log10(1 / mse_value)

def accuracy_smothness_error(g, f_hat, H, _lambda):
    data_accuracy = np.linalg.norm(g - ifft(fft(f_hat) * H))
    smothness = np.sum(np.square(np.gradient(f_hat)))
    return data_accuracy + _lambda * smothness

Wiener Filter:¶

$$\begin{align*} \hat f = \mathcal{F}^{-1} \Bigg\{\Bigg[\frac{H^*}{\|H\|^2 + K}\Bigg] G \Bigg\} \end{align*}$$

In [ ]:
def wiener_filter(G, H, K):
    W = np.conjugate(H) / (np.abs(np.square(H)) + K)
    F_hat = W * G
    f_hat = ifft(F_hat)
    return f_hat

Motion Blur Filter:¶

$$\begin{equation*} H(u,v) = T sinc\Big(\pi\big(u \Delta x + v \Delta y\big)\Big)\mathrm{e}^{\mathrm{i} \pi \big(u \Delta x + v \Delta y\big)} \end{equation*}$$

In [ ]:
def get_motion_blur_filter_fq(shape, T, dx, dy):
    U_DIM ,V_DIM = shape 
    H_u = np.arange(0, U_DIM, 1)
    H_v = np.arange(0, V_DIM, 1)
    
    H_U, H_V = np.meshgrid(H_v, H_u)
    H = np.dstack((H_U, H_V))

    dxy = np.array([dx, dy])

    def s(uv):
        return np.pi * uv.dot(dxy)
    
    H = np.apply_along_axis(s, axis=2, arr=H)
    H = T * np.sinc(H) * np.exp(-1j * H)    
      
    return H

Gausian Noise Filter:¶

$$\begin{align*} n(x,y) &\sim \mathcal{N}\big(\mu,\sigma^2\big)\\ N(u,v) &= \mathcal{F}\Big\{n(x,y)\Big\} \\ \end{align*}$$

In [ ]:
def get_gausian_noise_filter_fq(shape, variance, mean=0):
    n = np.random.normal(mean, variance, shape)
    N = fft(n)
    return N

Plotting:¶

In [ ]:
def plot_images_restauration(original_image, disturbed_image, restored_images, restored_titels):
    fig, ax = plt.subplots(2, math.ceil(len(restored_images) / 2) + 1)
    fig.set_size_inches((14, 12))

    ax[0,0].set_title(f'Original Image\n MSE: {mse(original_image, original_image):.2f}')
    ax[0,0].imshow(original_image, cmap='gray')
    ax[0,0].axis('off')

    ax[1,0].set_title(f'Disturbed Image\n MSE: {mse(original_image, disturbed_image):.2f}')
    ax[1,0].imshow(disturbed_image, cmap='gray')
    ax[1,0].axis('off')

    for i in range(len(restored_images)):
        x = int(i / 2) + 1
        y = i % 2
        ax[y, x].set_title(f'{restored_titels[i]}\n MSE: {mse(original_image, restored_images[i]):.2f}')
        ax[y, x].imshow(restored_images[i], cmap='gray')
        ax[y, x].axis('off')

    plt.tight_layout()
    plt.show()

Blind Wiener Filtering¶

Wiener Filtering With Optimization:¶

$$\begin{equation*} \underset{K}{argmin} \; \mathbb{E}\Big[\big\| \hat g - h \ast \hat f \big\|^2\Big] + \lambda \, \mathbb{E}\Big[\big\| \hat f^{\prime} \big\|^2\Big] \end{equation*}$$

In [ ]:
def optimize_k(disturbed_image_fq, motion_blur_filter_fq, _lambda=0, _dropout=True, _max_iter=50):

    disturbed_image = ifft(disturbed_image_fq)

    best_error = np.inf
    best_k = 0.5
    step = 0.5

    for _ in tqdm(range(_max_iter), desc='Optimizing K'):

        last_iter_error = best_error
        K = np.linspace(best_k + step, best_k - step, 5, endpoint=False)[1:]     
        
        for k in K:
            restored_image = wiener_filter(disturbed_image_fq, motion_blur_filter_fq, k)
            error = accuracy_smothness_error(disturbed_image, restored_image, motion_blur_filter_fq, _lambda)
            if error < best_error:
                best_error = error
                best_k = k

        step /= 2
        if _dropout and last_iter_error - best_error < 0.01:
            break

    return best_k

Wiener filtering By Averaging Over Noise:¶

[Yoo, Jae‐Chern, and Chang Wook Ahn. "Image restoration by blind‐Wiener filter." IET Image Processing 8.12 (2014): 815-823.]

Image-Restauration

In [ ]:
def wiener_filter_averaging(G, H, variance, n=10):
    f_hat = 0
    for _ in range(n):
        N_m = get_gausian_noise_filter_fq(G.shape, variance)
        D_hat_m = G - N_m
        S_D_hat_m= np.mean(np.abs(np.square(D_hat_m)))
        S_N_m = np.mean(np.abs(np.square(N_m)))

        f_hat_m = wiener_filter(D_hat_m, H, S_N_m/S_D_hat_m)
        f_hat += f_hat_m
    return f_hat / n
In [ ]:
def test_wiener_filter_averaging(image_path, _mb_T=1, _mb_dx=0.005, _mb_dy=0.001, _g_var=0.05):
    f = read_image(image_path)
    F = fft(f)
    H = get_motion_blur_filter_fq(F.shape, _mb_T, _mb_dx, _mb_dy)
    N = get_gausian_noise_filter_fq(F.shape, _g_var)
    G = H * F + N
    g = ifft(G)

    images = []
    titles = []
    variances = [0.01, 0.1, 0.2, 0.4, 0.6, 0.8]
    for var in tqdm(variances):
        images.append(wiener_filter_averaging(G, H, var))
        titles.append(f'Averaging (var): {var}')
    
    plot_images_restauration(f, g, images, titles)
In [ ]:
test_wiener_filter_averaging('images/schach.jpg', _mb_T=10, _mb_dx=0.01, _mb_dy=0.1, _g_var=0.1)

100%|██████████| 6/6 [00:13<00:00,  2.18s/it]
No description has been provided for this image

Other Posibilities Of Blind Wiener Filtering:¶

  • Parameter Estimation of Power Spectral Density of Noise and Image
  • location depentently K-parametric wiener filtering: K as adaptively varying-coefficient
  • ...

Testing:¶

In [ ]:
def test_restaurations(image_path, _mb_T=1, _mb_dx=0.005, _mb_dy=0.001, _n_var=0.05):
    f = read_image(image_path)
    F = fft(f)
    H = get_motion_blur_filter_fq(F.shape, _mb_T, _mb_dx, _mb_dy)
    N = get_gausian_noise_filter_fq(F.shape, _n_var)
    G = H * F + N
    g = ifft(G)

    # wiener filter with knowledge of powerspectrum of image and noise
    S_F = np.mean(np.abs(np.square(F)))
    S_N = np.mean(np.abs(np.square(N)))
    f_hat = wiener_filter(G, H, S_N/S_F)

    # wiener filter with optimization of K
    optimal_k = optimize_k(G, H, _lambda=0)
    f_hat_2 = wiener_filter(G, H, optimal_k)

    # wiener filter with averaging over noise
    f_hat_3 = wiener_filter_averaging(G, H, 0.01)

    f_hat_4 = wiener_filter_averaging(G, H, 0.2)

    plot_images_restauration(f, g, [f_hat, f_hat_2, f_hat_3, f_hat_4], ['Wiener Filter', 'Optimization', 'Averaging (var: 0.01)', 'Averaging (var: 0.2)'])

Only Motion Blur:¶

In [ ]:
test_restaurations('images/lion.jpg', _mb_T=1, _mb_dx=0.01, _mb_dy=0.01, _n_var=0)

Optimizing K:  56%|█████▌    | 28/50 [00:03<00:02,  8.81it/s]
No description has been provided for this image

Only Noise:¶

In [ ]:
test_restaurations('images/schach.jpg', _mb_T=1, _mb_dx=0, _mb_dy=0, _n_var=0.3)

Optimizing K:  26%|██▌       | 13/50 [00:17<00:48,  1.32s/it]
No description has been provided for this image

Motion Blur and Noise:¶

In [ ]:
test_restaurations('images/schach.jpg', _mb_T=1, _mb_dx=-0.001, _mb_dy=-0.02, _n_var=0.1)

Optimizing K:   6%|▌         | 3/50 [00:05<01:18,  1.67s/it]
No description has been provided for this image
In [ ]:
test_restaurations('images/eule.png', _mb_T=1, _mb_dx=0.001, _mb_dy=0.001, _n_var=0.2)

Optimizing K:   4%|▍         | 2/50 [00:07<02:56,  3.69s/it]
No description has been provided for this image
In [ ]:
test_restaurations('images/eule.png', _mb_T=10, _mb_dx=0.001, _mb_dy=-0.01, _n_var=0.01)

Optimizing K:  16%|█▌        | 8/50 [00:17<01:33,  2.23s/it]
No description has been provided for this image
In [ ]:
test_restaurations('images/schach.jpg', _mb_T=10, _mb_dx=0.1, _mb_dy=0.2, _n_var=0.1)

Optimizing K:   6%|▌         | 3/50 [00:04<01:14,  1.59s/it]
No description has been provided for this image